library(fixest)
library(modelsummary)
library(spmoran)
library(REndo)
library(WeightIt)
library(DirectEffects)
library(vtable)
library(spdep)
library(panelView)
library(ggplot2)
library(ggpubr)
library(viridis)
library(patchwork)
library(cragg)
library(mediation)
library(sensemakr)
library(robomit)
library(spatialreg)
library(rdrobust)


data <- readRDS("data_final_merged.rds")  %>% dplyr::select(NUTS_ID, Judet, Year, Month, p_ottoman, share, count, ps_1699, lon, lat, tavg, lograin, logelev, logslope, wheat_rainfed, log_large_rivers_distance,log_sea_lines_distance,route_density_2, vienna_distance, istanbul_distance,geometry)

data <- data %>%
  mutate(after_feb_2020 = ifelse(Year > 2020 | (Year == 2020 & Month > 2), 1, 0))

model1 <- feols(share ~ p_ottoman* after_feb_2020 , data = data, vcov = "conley")
model2 <- feols(share ~ p_ottoman* after_feb_2020 + count , data = data, vcov = "conley")
model3 <- feols(share ~ p_ottoman* after_feb_2020 + count| Year + Month , data = data, vcov = "conley")
model4 <- feols(share ~ p_ottoman* after_feb_2020 + count + lat + lon | Year + Month , data = data, vcov = "conley")
model5 <- feols(share ~ p_ottoman* after_feb_2020 + count + tavg + logelev + logslope + lograin + log_large_rivers_distance + log_sea_lines_distance| Year + Month,  data = data, vcov = "conley")

msummary(list(model1, model2, model3,model4,model5), stars = c('.' = .1,'*' = .5, '**' = .01, '***' = .001))

